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We determine the first independent part of tfie coefficient in tfie weak coupiing expansion of 
the QCD pressure at high temperatures, the one proportionai to the maximai power of the number 
of quark fiavors A'f. In addition to introducing and deveioping computationai methods that can 
be used in evaluating other parts of the expansion, our calculation provides a result that becomes 
dominant in the limit of large A'f and a fixed effective coupling gig = g'^Nf/2. 



PACS numbers: ll.lO.Wx, 11.15.1 



I. INTRODUCTION 
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The single most fundamental quantity characterizing 
the bulk properties of a system in thermodynamic equi- 
librium is its partition function, or equivalently the func- 
tional dependence of its pressure p on the temperature T 
(and other parameters such as the chemical potentials). 
In the case of hot quantum chromodynamics (QCD), this 
function has been extensively studied both on the lattice 
[U, 01 and using perturbation theory, with the purpose of 
obtaining a consistent description of the quantity all the 
way from the phase transition region [T ~ Tc) to asymp- 
totically high temperatures. At present, the state of the 
art on the perturbative (high T) side is order In g in 
the strong coupling constant, which has been reached 
both at zero and finite quark chemical potentials 0, [J ■ 
In addition, perturbative methods have been successfully 
applied to a plethora of other quantities (for some recent 
results, see e.g. Rcfs. [5, 6, 7] and references therein). 

A persistent problem in perturbative finite tempera- 
ture QCD is the slow convergence of the various weak 
coupling expansions. For most quantities, the issue can 
be traced back to the contributions of the soft {gT) and 
ultrasoft {g^T) energy scales, which respectively enter 
through the three-dimensional effective theories EQCD 
and MQCD Q. In the case of the pressure, the contribu- 
tions of these scales to the term in the expansion 
have by now been determined [1, HH, , and to this 
end, there is clear motivation for evaluating the part cor- 
responding to the hard (27rT) scales as well This 
requires the evaluation of all four-loop vacuum diagrams 
in the full four-dimensional theory, a task that has al- 
ready been completed in the similar, thou gh t echnically 
much simpler, case of scalar theory [l^ Il5l|. 

The 0{g^) contribution of the QCD pressure can be di- 
vided into several gauge invariant parts, proportional to 
various group theory invariants, most importantly pow- 
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ers of the number of fundamental fermion flavors TVf . In 
the present paper, our purpose is to evaluate the first — 
and computationally most straightforward — of these, 
the one proportional to the maximal power Nf . Our 
motivation for this is twofold. On one hand, we wish 
to demonstrate that the computational machinery built 
for three-loo p Q CD and four-loop 0'* theory calculations 

m Refs. [11, Mini 

can be straightforwardly applied to 
the four-loop level of QCD as well. In addition, the Nf 
term not only represents the first independent piece of the 
0{g^) coefficient, but in fact becomes dominant over all 
the other contributions in the limit of a large flavor num- 
ber, where the effective coupling = Tpg'^ = Nfg'^/2 is 
kept fixed while N{ is taken to infinity. In this limit, the 
theory in fact somewhat trivializes, enabling an all orders 
numerical evaluation of the partition function [l^: 
and providing a rough numerical check for our result. 

The paper is organized as follows. In the rest of the 
first Section, we present our notation and explain, how 
renormalization is performed in our work. After this, we 
proceed to review the organization of our calculation in 
Section II, while the bulk of the detailed computations 
is left to Section III. Our result is finally assembled and 
discussed in Section IV, where we in addition analyze 
the convergence of the weak coupling expansion of the 
pressure and draw our final conclusions. 



A. Notation 

We work in 0? = 4 — 2e dimensional Euclidean space- 
time, using dimensional regularization to regulate both 
ultraviolet (UV) and infrared (IR) divergences. We use 
lower case bold letters to denote three-dimensional spa- 
tial vectors and capital letters for four-dimensional space- 
time vectors P ~ {po,p), so that 
As usual, sum-integrals are defined by 



PiPi 



A' 



P 



T 



i2n) 



d-l 



An 



{2n: 



(1) 



' ^tJ2 h (2) 



2 



where A is the MS renormahzation scale and A the MS 
one, and P/{P} refer to bosonic/fermionic Matsubara 
modes, respectively. Bosonic sum-integrals, in which the 
Po = term has been subtracted out, are marked with a 
prime, and, borrowing notation from Refs. [8,, Jii], some 
standard sum-integrals are denoted by 



(P2)«' 
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We perform renormahzation by consistently using the 
bare coupling f/s in all of our calculations, and only in 
the end expressing the result in terms of the physical, 
renormalized one via the relation = ZgCj^. We need 
the parameter Zg to two-loop order, to which it reads 
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Here, the coefficients Pi can for our purposes be deduced 
from the known exact beta- function of A'f — > oo QCD, 
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In this limit, the running coupling correspondingly reads 
1 1 Tpg^fi), A' 



1 



67r2 
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which this time is exact to all loop orders. It exhibits the 
well-known Landau singularity, which implies that large- 
Nf QCD is well-defined only below some energy scale 
Al- 

Finally, in the course of our calculation, we will often 
switch to dimensionless integration variables without ex- 
plicitly saying so. In all of these cases, the integration 
momenta and coordinates have been scaled by the ap- 
propriate power of 27rr. 



II. ORGANIZING THE CALCULATION 

The weak coupling expansion of the QCD pressure is known to order In g for an arbitrary number of colors and 
massless fundamental fermion flavors 0,11]. In the Nf oo limit (with fixed ^2^^) and at zero chemical potential, the 
result reads [2ll 
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in which we wish to compute the unknown coefficient 
pg that is nothing but the Njg^ term in the full theory 
pressure. Here, we have subtracted out the pressure of 
a gas of free quarks and gluons, po — ir^T'^/Ab {(Ia + 
7NcNf/4), where (Ia = — 1 and which dominates the 
pressure in the large N{ limit (for all .g^ff). 

At leading order in the large N{ expansion oip—po, the 
effective three-dimensional theories EQCD and MQCD, 
which incorporate the effects of the soft scales gT and 
f;2p to the pressure, do not contribute to the coefficients 
of even powers of g. We can thus concentrate on evaluat- 
ing the contribution of the hard scale, corresponding to 
the strict (unresummed) perturbation expansion of the 
pressure in the four-dimensional theory, commonly de- 
noted by pe- We note that while this function suffers 
from IR divergences that for finite N( are canceled by 
UV divergences in the effective theory contributions, in 
the large A^f limit they are all guaranteed to vanish in 



dimensional regularization. 

In the large N{ limit, the function pE has a simple 
diagrammatic expansion of the form 



dA 2^435^ 



where we have denoted 
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and the integrals correspond to the three Feynman 
graphs in Fig. [T] Here, the one-loop large iVf gluon self 
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FIG. 1: The two-, three- and four-loop Feynman diagrams contributing to pE in the large A'f limit. The straight lines with an 
arrow correspond to quarks, and the wavy lines to gluons. 



energy reads 

(2Q_P)^{2Q ~ P), 

'{Q} 



(16) 



in which the color indices have been suppressed (the trace 
over them having already been carried out). 

At this point, we may simplify the calculations by di- 
viding the self energy into its three-dimensionally lon- 
gitudinal and transverse components (for a precise def- 
inition, see e.g. Ref. (2^). in terms of which the above 
integrals read 
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Extracting now the UV divergences of the self-energies 
by denoting 



IiL,T{P) = ^L/t{P) 
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the expression for the function pE up to order obtains 
the remarkably simple form 
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Here, the bars in /f,.c denote the replacing of I{]^/t{P) 
by ni/7-(P) (note that with this definition la = la), 
and Eq. ([9|) has been used to cancel terms originating 
from the renormalization corrections against those left 
over from the redefinition of Eq. ([^0]) . As la and /& are 
both known, we are only left with the task of computing 
Ic, which is furthermore observed to be finite in dimen- 
sional regularization. 

For notational simplicity, it is convenient to separate 
the r = contributions of the self-energies by writing 
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We find that, in various limits, the two parts of the func- 
tions can be written as: 
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where (n,p) = (po/(27rT), p/(27rT)), 

• at P 00 (up to but not including order / P^) 
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at p ^ (up to but not including order p'^): 
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The required integral can then be divided into four 
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parts, 
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The values of all of the integrals appearing here are avail- 
able in the literature (see e.g. Ref. [16]), giving in the end 
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Of these terms, Ki and are UV divergent but IR 
finite, while and K4 are UV finite but IR divergent. 
In the following Section, we will set out to evaluate these 
functions one by one. 



III. EVALUATION OF THE INTEGRALS 

In performing the integrals K1-K4, we employ the gen- 
eral strategy of first identifying their UV and IR diver- 
gent parts, and then separating them from the rest, writ- 
ing 



^ + iff". 



(31) 



The divergent pieces must be evaluated analytically in 
rf = 4 — 2e dimensions, whereas in the finite integrals we 
may set d — -i and use the integral representations of 
Eqs. ((24|) and (|25|) for the self-energies. Apart from the 
case of K4, which is treated separately, the subtracted 
divergent pieces are always identified with specific parts 
of the integrands of Eqs. (|24|) and (|25|) . corresponding to 
terms in the Taylor expansion of the hyperbolic cosecant. 

For ease of notation, we write the finite parts of the 
Ki in terms of dimensionless integrals Mi, 
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which one can, when necessary, easily evaluate with nu- 
merical tools. 



A. Ki integral 

Using the expression of Eq. (^5]) for I^'^\p), the in- 
tegral Ki can be trivially expressed in terms of the X™ 
functions. 
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B. K2 integral 

The K2 integral is logarithmically UV divergent due to 
the integrand behaving like 1/ at large P. Separating 
thus its large P limit from the rest, we easily find for the 
divergent part 
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where all terms are again readily available. 

The finite part of the integral is most conveniently eval- 
uated in two pieces. Starting from the Po 7^ terms, we 
obtain using the notation of Eq. ([5^ 



M' 



n=l 



d^p /5 , A^ 2 ix 



e icscnr ^ + g 



5 , A2 
— h m — 7: — 
3 47r2T2 



.,'5 , A2 
-2(-+ln 



47r2T2 ' 2 



/2 
2 1 



(37) 



where we have defined 
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Here, the last term inside the parentheses is a result of 
the UV subtraction we have performed. 

Concentrating on the three integrals to be evaluated, 
we immediately observe that = 0, as the p integral 
in this case yields S{r), while the rest of the integrand 
vanishes at r = 0. For M2^, we on the other hand note 
that using 
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both the p integral and the Matsubara sum can be easily 
performed, leaving us with a straightforwardly solvable 
hyperbolic integral with the result 
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Finally, the remaining integral M2 we evaluate numeri- 
cally, obtaining 
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0.0099981763. 
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Collecting all the different parts of our result, we finally 
find for K2 
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Contrary to the po 7^ case, the pa — Q contribution 
to K2 is UV finite, and the UV subtraction term in fact 
yields zero in dimensional regularization. The finite di- 
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we may write Eq. in the form 
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We then divide the integral into three parts in a manner 
analogous to Eq. ([57)) . 
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The evaluation of these integrals is easily completed. 
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C. integral 

The K^g, integral is UV finite, but its po — Q term con- 
tains an IR divergence and must therefore be considered 
separately. There, the divergent part is easily isolated 
using Eq. ((28| and is given by 
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This integral, however, obviously vanishes in dimensional 
regularization. 

In the IR finite part of the pq — term, we can set 
e = and thus get 
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where the last term corresponds to the subtraction of the 
divergence. This yields 
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With the IR finite po ^ part, we write analogously 
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For Mg", we get after integrating over p and summing 
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while the remaining integral, Afg^, must again be evalu- 
ated numerically with the result 
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Summing up all the different pieces, we finally get for 
the entire integral 
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With the integral K^, we again have an IR divergence 
in the = term, which we attempt to separate first. 
Using Eq. ([Ml), we get 
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which, just like K^^^ , vanishes in dimensional regulariza- 
tion. 

For the IR finite part of the po — term, we on the 
other hand obtain after subtracting off all the IR diver- 
gent parts 
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Here, we identify the potential problem that as the two 
first terms inside the square parentheses are separately 
divergent and only converge when summed together, the 
numerical evaluation of the integral is inconvenient in its 
present form. To this end, we consider the first term in 
more detail. 

By carrying out the angular integrals, we find for the 
function i?+(0,p) 
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to which we add and subtract a (p-independent) term in 
such a way that its small p behavior becomes apparent. 
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where the first term obviously cancels against the second 
term of Eq. ((64| . 

Putting everything together, we obtain 
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In all of these formulas, the r, s and t derivatives are understood to act only on the sin{px)/{px) terms. 



Finally, the po ^ piece of K4, M^, is IR finite and in fact closely analogous to M, 
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FIG. 2: [color online] The behavior of the different orders of the weak coupling expansion of the large A'f pressure, Eq. 
normalized to the pressure of free gluons po,giuon = dA7r^r*/45 and plotted together with the all orders numerical result of 
Ref. [l^. All the results are evaluated with the renormalization scale A = e~'^^nT, corresponding to the canonical 'optimal' 
scale of Ref. [2^. For the curve, we also display the effects of varying the renormalization scale by a factor e"^^ « 1.78, which 
corresponds to the shaded grey area. Note that when viewing the pressure as a function of the coupling, also the all orders 
numerical result has an explicit dependence on the renormalization scale, which only vanishes when taking the running of g^ft 
into account. This implies that the grey area is not to be compared to the all orders result, but is displayed only to give an 
idea of the size of the scale dependence of the perturbative result. 



Pulling all the pieces together, we find for the entire K4 integral 



K4 



135(27r)4 



X 27.55287622. 



(72) 



r 



IV. DISCUSSION AND CONCLUSIONS 

In the previous Section, we have evaluated all of the 
integrals Ki needed in constructing a result for the func- 
tion pq of Eq. (fTTjl . Collecting everything with the help of 
Eqs. ([2T|) and (|30p , we now finally obtain for the quantity 



800, 2 A 
In 

9 AttT 



16 

y 



80C'(-1) 40C'(-3) 

y c(-i) ~ y c(-3) 



f 207E 

ln> 



;ln2 



(73) 



A-kT 



296.055373. 



From here, we may also read off the leading large A^f 
behavior of the function qdNf) defined in Eq. (6.11) of 
Ref. i, 

q^Nf) = 0.469861169 + 0{N^), (74) 
as well as the coefficient Cg defined in Ref. [20l |. 

Ce = 21.8597689. (75) 



It is quite reassuring that our result turned out to lie 
within the error bars of the very non-trivial numerical 
estimate Cg = 20(2) of Ref. 0. 

The above result can be contrasted on one hand with 
the exact numerical solution of the large A^f pressure [HI , 
and on the other hand to the previous terms in the expan- 
sion of Eq. (jlip . In Fig. [21 we perform both comparisons, 
displaying the perturbative result for the pressure to var- 
ious orders as a function of the effective coupling gcs, and 
comparing it at the same time to the all orders resum- 
mation of Ref. [l^ . As noted already in previous studies 
[19L [20I I , the convergence of the weak coupling expansion 
seems quite impressive: The region of applicability of the 
result clearly increases order by order, and it is only for 
relatively large values of the effective coupling, g^g > 20, 
that the renormalization scale dependence of the order 
result becomes strong enough to completely ruin its 
predictive power. 

The good convergence properties of the large A^f pres- 
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sure seem to indicate that in this case, the smah parame- 
ter in the weak couphng expansion is really g^ff/(47r)^, as 
can indeed be verified from Eq. pT|) . This should be con- 
trasted to the case of full QCD, where the EQCD contri- 
butions come with an associated expansion parameter g, 
while the MQCD ones are completely non-perturbative. 
Remarkably, in the large Nf limit, even the effective the- 
ory contributions organize themselves in the form of g^f^ 
times a series expansion in f;^fj/(47r)^. This can be seen 
to follow from the fact that in the large Ni limit, EQCD 
is a free theory [23|. Its contribution to the pressure has 
the form m| x /((/l/mE), where ttie ~ 5 is the mass 
parameter and ~ g"^ the gauge coupling constant of 
EQCD, both having weak coupling expansions in pow- 
ers of g^/(47r)^. In the large Nf limit, the ratio gl/ruE 
behaves as 1/N{, implying that the pressure of EQCD 
is then indeed merely a number times rrig, and there- 
fore contributes to the pressure of the full theory as g^^f 
times an expansion in g'^^ / (Att)'^ . Finally, the pressure of 
MQCD is entirely subleading in the large N{ limit. 

Apart from providing the coefficient of the N^g^ 
term in the expansion of the QCD pressure, our work 
has demonstrated the applicability of the computa- 
tional methods developed for lower order calculations in 
Refs. [TgI . [iTj and for the case of scalar <j)^ theory in 
Ref. [ij] in tackling four- loop computations in QCD. We 
have been able to perform our work semi-analytically, 
dealing with all UV and IR divergent parts and loga- 
rithms of the renormalization scale analytically in dimen- 



sional regularization, and evaluating the remaining finite 
parts numerically to a high accuracy. At worst, we have 
had to deal with two-fold numerical integrals, which is 
a very modest challenge for state of the art computing 
facilities. 

The success of the calculation performed in this paper 
encourages one to pursue further pieces of the coeffi- 
cient of the QCD pressure using similar machinery. In 
particular, the term proportional to Nfg^ in the result 
consists of a reasonably limited set of diagrams, a large 
number of which seem amenable to an analogous treat- 
ment. Beyond this, the task rapidly becomes much more 
complicated, and an approach based on computer algebra 
methods proves necessary. 
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